      program btu0
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      parameter (np1=4,mp1=np1+1)
      parameter (ftol=1.d-10,nrep=999)
c
      external famoeb,gamoeb
      character*6 name(np1)
      character*24 fdate
      dimension x1(np1),p1(mp1,np1),y1(mp1),x0(np1)
      dimension is(nplay,ndec)
      dimension est1(np1),xx1(np1)
      dimension bllf1(nrep),boot1(np1,nrep),bratio(nrep),
     &          btmp(nrep),bs1(np1),bss1(np1)
      common /mats/is
      common /est/est1
c
      data name/'BETA','PZ','SB','SZ'/
c
      open(unit=11,file='ra2011.data')   !with card sleeves on boards
      open(unit=10,file='btu0.dat')
c
      read(11,*)((is(i,n),n=1,ndec),i=1,nplay)
c
c-----------------------------------------------------------------
      write(*,'('' ***** btu0 ******'',20x,a24,
     &//"    Risk & Ambiquity Experiment")')fdate()
      write(*,'(/"  All 2011 Sessions"//)')
c
      call fgrid
      call fch
c
c Set parameters
         x0(1) = dlog(1.d0/0.9d0 - 1.d0)
         x0(2) = dlog(0.5d0/0.2d0 - 1.d0) 
         x0(3) = dlog(1.d0/0.6d0 - 1.d0)
         x0(4) = dlog(1.d0/0.15d0 - 1.d0)
c
         delta=0.1d0
         iter=1000  
         call init(x0,p1,y1,famoeb,delta)
         call amoeba(p1,y1,mp1,np1,np1,ftol,famoeb,iter)
         call avg(p1,x1,np1,mp1)
         z1 = -famoeb(x1)
         call m2h(x1,est1,np1)
         write(*,'("est = ",4g13.5)')(est1(k),k=1,np1)
         write(*,'(/"Initial LL = ",g14.6)')z1
         write(10,910)0,0.0,z1,(est1(k),k=1,np1)
c         stop
c 
         call savephi(est1)
c
c begin interations
c
      dseed = 90001.d0
      do 500 it=1,nrep
   25 xseed = dseed
         call data(it,dseed)
         call fch
c
c uni model

         delta=0.1d0
         iter=2000
         call init(x0,p1,y1,famoeb,delta)
         call amoeba(p1,y1,mp1,np1,np1,ftol,famoeb,iter)
         if (iter.ge.1000)  then
           write(*,'(5x,"irep =",2x,i5)')it
           go to 25
         endif
         call avg(p1,x1,np1,mp1)
         bllf1(it) = -famoeb(x1)
         call m2h(x1,xx1,np1)
         do 30 k=1,np1
            boot1(k,it) = xx1(k)
   30    continue
c
         write(10,910)it,xseed,bllf1(it),(boot1(k,it),k=1,np1)
  500 continue
  910   format(i4,2x,f20.1,5(2x,g14.6))
  911   format(2x,10(2x,g14.6))
c
c --- process the bootstrap results ---
c
      do 160 j=1,np1
         bs1(j)=0.d0
         bss1(j)=0.d0
  160 continue
         blm=0.0d0
         blms=0.0d0
      do 162 k=1,nrep
         blm=blm+bllf1(k)
         blms=blms+bllf1(k)*bllf1(k)
         do 162 j=1,np1
            bs1(j)=bs1(j)+boot1(j,k)
            bss1(j)=bss1(j)+boot1(j,k)*boot1(j,k)
  162 continue
      blm=blm/float(nrep)
      blms=(blms-float(nrep)*blm*blm)/float(nrep-1)
      blms=dsqrt(blms)
      do 165 j=1,np1
         bs1(j)=bs1(j)/float(nrep)
         bss1(j)=(bss1(j)-float(nrep)*bs1(j)*bs1(j))/float(nrep-1)
         bss1(j)=dsqrt(bss1(j))
  165 continue
      write(*,810)nrep
      t1 = 0.05*float(nrep)
      nlo = t1 + 1
      nhi = nrep - nlo + 1
      do 167 j=1,np1
         tratio=est1(j)/bss1(j)
         do 166 k=1,nrep
  166       btmp(k)=boot1(j,k)
         call sort(nrep,btmp)
         xlo=btmp(nlo)
         xhi=btmp(nhi)
         write(*,811)name(j),bs1(j),bss1(j),tratio,xlo,xhi
         write(*,813)est1(j)
  167 continue
      write(*,'(//"The mean and std dev of log-likelihood:")')
      write(*,812)blm,blms
  810 format(//'*** BOOTSTRAP RESULTS ***  nrep=',i5//
     &  'name ','  estimate  ','  std dev   ','  t-ratio   ',
     &      2x,'---95 percent conf. int.--- ')
  811 format(a6,2x,2g13.6,2x,2g11.4,x,2g13.6)
  812 format(/2g14.6)
  813 format(7x,g14.6)
c
      end
c***********************************************************************
      subroutine m2h(x,xx,np)
      implicit real*8 (a-h,o-z)
      dimension x(np),xx(np)
c
      xx(1) = 1.d0/(1.d0+dexp(x(1)))
      xx(2) = 0.5d0/(1.d0+dexp(x(2)))+0.5d0
      xx(3) = 1.d0/(1.d0+dexp(x(3)))+1.d-20
      xx(4) = 1.d0/(1.d0+dexp(x(4)))+1.d-20
c
      return
      end
c***********************************************************************
      double precision function famoeb(x)
      implicit real*8 (a-h,o-z)
      parameter (np=4)
      dimension x(np)
c
      beta = 1.d0/(1.d0+dexp(x(1)))
      pz = 0.5d0/(1.d0+dexp(x(2)))+0.5d0
      sb=1.d0/(1.d0+dexp(x(3)))+1.d-5
      sz=1.d0/(1.d0+dexp(x(4)))+1.d-5
c
      call fam(beta,pz,sb,sz,pt)

      famoeb = -pt
c
      return
      end
c***********************************************************************
      subroutine fam(beta,pz,sb,sz,pt)
      implicit real*8 (a-h,l,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pch(nplay,ngrid,ngrid),phi(ngrid,ngrid)
      common /pchm/pch

      call fprob(beta,pz,sb,sz,phi)
c
      pt = 0.d0
      do 50 id=1,nplay
         pid = 0.d0
         do 20 j=1,ngrid
            do 15 k=1,ngrid
               pid = pid + phi(j,k)*pch(id,j,k)
   15       continue
   20    continue
         pid = dlog(pid)
         pt = pt + pid
   50 continue
      return
      end
c***********************************************************************
      subroutine fgrid
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension w0(ngrid),w(ngrid,ngrid)
      dimension px(ndec,2),pxm(ngrid,ngrid,ndec,2)
      common /grid/pxm,w
c
      do 5 i=1,ngrid
         w0(i) = 1.d0
         if((i.eq.1).or.(i.eq.ngrid)) w0(i) = 0.5d0
    5 continue
c
      dz=1.d0/float(ngrid-1)
      do 30 j=1,ngrid
         beta = dz*(j-1)
         do 20 k=1,ngrid
            if (k.lt.ngrid) then
               pz = 0.5d0 + 0.5d0*dz*(k-1)
            else
               pz = 0.999d0
            endif
            gam = 2.d0*dlog(pz/(1.d0-pz))
            call fpx(gam,beta,px)
            do 10 m=1,ndec
               pxm(j,k,m,1) = px(m,1)
               pxm(j,k,m,2) = px(m,2)
   10       continue
            w(j,k) = w0(j)*w0(k)
   20    continue
   30 continue
c
      return
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fch
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension pch(nplay,ngrid,ngrid)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /mats/is
      common /pchm/pch

      do 50 id=1,nplay
         do 30 j=1,ngrid
            do 25 k=1,ngrid
               pin = 1.d0
               do 20 m=1,ndec
                  ix = is(id,m)
                  pin = pin*pxm(j,k,m,ix)
   20          continue
               pch(id,j,k) = pin
   25       continue
   30    continue
   50 continue
      return
      end
c***********************************************************************
      subroutine fprob(beta,pz,sb,sz,phi)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,ngrid=51)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension phi(ngrid,ngrid)
      common /grid/pxm,w
c
      pt = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 20 j=1,ngrid
         bj = dz*(j-1)
         pb = (bj - beta)/sb
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz)/sz
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            phi(j,k) = w(j,k)*pg*pb
            pt = pt + phi(j,k)
   10    continue
   20 continue
      do 30 j=1,ngrid
         do 25 k=1,ngrid
            phi(j,k) = phi(j,k)/pt
   25    continue
c      write(*,'(11f7.3)')(phi(j,k),k=1,ngrid)
   30 continue
c      stop
      return
      end
c***********************************************************************
      subroutine savephi(x)
      implicit real*8 (a-h,o-z)
      parameter (np=4,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension x(np),gl(npts)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      common /grid/pxm,w
      common /glm/gl
c
      beta = x(1)
      pz = x(2)
      sb = x(3)
      sz = x(4)
      write(*,'(4g12.4)')beta,pz,sb,sz
c
      gt = 0.d0
      dz = 1.d0/float(ngrid-1)
      do 20 j=1,ngrid
         j0 = ngrid*(j-1)
         bj = dz*(j-1)
         pb = (bj - beta)/sb
         pb = pb**2
         pb = dexp(-0.5d0*pb)
         do 10 k=1,ngrid
            pk = 0.5d0 + 0.5d0*dz*(k-1)
            pg = (pk - pz)/sz
            pg = pg**2
            pg = dexp(-0.5d0*pg)
            gl(j0+k) = w(j,k)*pg*pb
            gt = gt + gl(j0+k)
   10    continue
   20 continue
c
      gl(1) = gl(1)/gt
      do 30 i=2,npts
         gl(i) = gl(i-1) + gl(i)/gt
   30 continue
c      write(*,'(21f8.4)')(gl(i),i=1,npts)
c      stop
      return
      end
c***************************************************************
      subroutine data(it,dseed)
      implicit real*8 (a-h,o-z)
      parameter (nplay=101,ndec=12,ngrid=51,npts=ngrid*ngrid)
      dimension pxm(ngrid,ngrid,ndec,2),w(ngrid,ngrid)
      dimension gl(npts)
      dimension is(nplay,ndec)
      common /grid/pxm,w
      common /glm/gl
      common /mats/is
c
      do 50 id=1,nplay
         r = ggubfs(dseed)
         do 15 i=1,npts
            if (gl(i) .gt. r) goto 16
   15    continue
c
   16    j = (i-1)/ngrid + 1
         k = i - ngrid*(j-1)
c
         do 30 m=1,ndec
            x = ggubfs(dseed)
            p1 = pxm(j,k,m,1)
            if (x .le. p1) then 
               is(id,m)=1
            else
               is(id,m)=2
            endif
   30    continue
c      write(*,'(13i3)')id,(is(id,k),k=1,ndec)
   50 continue
      return
      end
c
c***********************************************************************
      subroutine init(xint,p,y,famoeb,delta)
      implicit real*8 (a-h,o-z) 
      parameter (np=4,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-famoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
c      delta=0.1d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end
c**************************************************************
      subroutine sort(n,ra)
      implicit real*8 (a-h,o-z)
      dimension ra(n)
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
        l=l-1
        rra=ra(l)
        else
        rra=ra(ir)
        ra(ir)=ra(1)
        ir=ir-1
        if(ir.eq.1)then
        ra(1)=rra
        return
        endif
        endif
        i=l  
        j=l+l
20      if(j.le.ir)then
        if(j.lt.ir)then
        if(ra(j).lt.ra(j+1))j=j+1
        endif
        if(rra.lt.ra(j))then
        ra(i)=ra(j)
        i=j
        j=j+j
        else
        j=ir+1
        endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end      
c========================================================================
      double precision function ggubfs (dseed)
      double precision   dseed
      double precision   d2p31m,d2p31
      data              d2p31m/2147483647.d0/
      data              d2p31 /2147483648.d0/
      dseed = dmod(16807.d0*dseed,d2p31m)
      ggubfs = dseed / d2p31
      return
      end

